// Copyright 2010 Google Inc.
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//      http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
//
// Author: tu@tulrich.com (Thatcher Ulrich)

// C++ translation of the WebGL benchmark code written in Javascript.
// Math utilities.

#include <math.h>
#include "glmath.h"

void CrossProduct(const float a[3], const float b[3], float out[3]) {
  // a1b2 - a2b1, a2b0 - a0b2, a0b1 - a1b0
  out[0] = a[1] * b[2] - a[2] * b[1];
  out[1] = a[2] * b[0] - a[0] * b[2];
  out[2] = a[0] * b[1] - a[1] * b[0];
}

float DotProduct(const float a[3], const float b[3]) {
  return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}

void VectorNormalize(const float v[3], float out[3]) {
  float l2 = DotProduct(v, v);
  if (l2 <= 0) {
    // Punt.
    out[0] = 1;
    out[1] = 0;
    out[2] = 0;
    return;
  }
  float scale = 1 / sqrtf(l2);
  out[0] = v[0] * scale;
  out[1] = v[1] * scale;
  out[2] = v[2] * scale;
}

void SetAffineMatrix(float e0, float e4, float e8, float e12,
                  float e1, float e5, float e9, float e13,
                  float e2, float e6, float e10, float e14, float out[16]) {
  out[0]  = e0;
  out[4]  = e4;
  out[8]  = e8;
  out[12] = e12;
  out[1]  = e1;
  out[5]  = e5;
  out[9]  = e9;
  out[13] = e13;
  out[2]  = e2;
  out[6]  = e6;
  out[10] = e10;
  out[14] = e14;

  out[3] = 0;
  out[7] = 0;
  out[11] = 0;
  out[15] = 1;
}

void SetMatrix(float e0, float e1, float e2, float e3,
               float e4, float e5, float e6, float e7,
               float e8, float e9, float e10, float e11,
               float e12, float e13, float e14, float e15,
               float out[16]) {
  out[0] = e0;
  out[1] = e1;
  out[2] = e2;
  out[3] = e3;

  out[4] = e4;
  out[5] = e5;
  out[6] = e6;
  out[7] = e7;

  out[8] = e8;
  out[9] = e9;
  out[10] = e10;
  out[11] = e11;

  out[12] = e12;
  out[13] = e13;
  out[14] = e14;
  out[15] = e15;
}

void MakeOrientationAffine(const float pos[3], const float dir[3],
			   const float up[3], float model_mat[16]) {
  float right[3];
  CrossProduct(dir, up, right);

  SetAffineMatrix(
      dir[0], up[0], right[0], pos[0],
      dir[1], up[1], right[1], pos[1],
      dir[2], up[2], right[2], pos[2],
      model_mat);
}

// As above, but puts result in given output object.
void TransformPoint(const float t[16], const float p[3], float out[3]) {
  out[0] = t[0] * p[0] + t[4] * p[1] + t[8]  * p[2] + t[12];
  out[1] = t[1] * p[0] + t[5] * p[1] + t[9]  * p[2] + t[13];
  out[2] = t[2] * p[0] + t[6] * p[1] + t[10] * p[2] + t[14];
}

// Return the inverse of the rotation part of matrix a, assuming
// that a is normalized.  This is just the transpose of the 3x3
// rotation part.
void InvertNormalizedRotation(const float a[16], float out[16]) {
  SetAffineMatrix(
      a[0], a[1], a[2], 0,
      a[4], a[5], a[6], 0,
      a[8], a[9], a[10], 0,
      out);
}

// Return the inverse of the given affine matrix, assuming that
// the rotation part is normalized, by exploiting
// transpose==inverse for rotation matrix.
void InvertNormalized(const float m[16], float out[16]) {
  InvertNormalizedRotation(m, out);
  float trans_prime[3];
  float p[3] = { m[12], m[13], m[14] };
  TransformPoint(out, p, trans_prime);
  out[12] = -trans_prime[0];
  out[13] = -trans_prime[1];
  out[14] = -trans_prime[2];
}

void MakeViewFromOrientation(const float orient[16],
			     float view[16]) {
  float temp[16];
  SetAffineMatrix(
      orient[8],  orient[4], -orient[0], orient[12],
      orient[9],  orient[5], -orient[1], orient[13],
      orient[10], orient[6], -orient[2], orient[14],
      temp);
  InvertNormalized(temp, view);
}

// fovy in radians
// aspect ratio is x/y
void MyGluPerspective(float fovy, float aspect, float z_near, float z_far,
                      float out[16]) {
  float f = 1 / tanf(fovy / 2);
  SetMatrix(
      f / aspect, 0, 0, 0,
      0, f, 0, 0,
      0, 0, (z_far + z_near) / (z_near - z_far), -1,
      0, 0, 2 * z_far * z_near / (z_near - z_far), 0,
      out);
}
 
void MyMatrixMultiply(const float ma[16], const float mb[16], float out[16]) {
  for (int i = 0; i < 4; i++) {
    for (int j = 0; j < 4; j++) {
      float val = 0;
      for (int k = 0; k < 4; k++) {
	val += ma[j + k * 4] * mb[k + i * 4];
      }
      out[j + i * 4] = val;
    }
  }
}

